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experiment,  in  which  the  movement  of  cabbage  root  flies  in  the  presence  of  a 
cabbage  crop  vas  studied.  It  is  proposed  that  such  a  parameter  estimation 
method  can  be  a  useful  analytical  tool  to  help  develop  appropriate  models  in 
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SPLINE-BASED  PARAMETER  ESTIMATION 
TECHNIQUES  FOR  TVO-DIMENSIONAL 
CONVECTION  AND  DIFFUSION  EQUATIONS 

L. L.  Zia 

ABSTRACT 

A  general  approximation  framework  based  on  bicubic  splines  is 
developed  for  estimating  temporally  and  spatially  varying  parameters  in 
two-dimensional  convection  and  diffusion  equations  derived  from  mass  transport 
theory.  The  "parameter  estimation  problem"  is  first  cast  as  an  abstract 
infinite  dimensional  minimization  problem.  Then  a  sequence  of  approximate, 
finite  dimensional  minimization  problems  is  defined,  which  yields  a  sequence 
of  parameter  estimates.  Convergence  results  relating  the  approximate  problems 
to  the  full  infinite  dimensional  problem  are  presented,  as  well  as  a 
discussion  addressing  computer  implementation.  Finally,  the  technique  is 
applied  to  the  analysis  of  actual  biological  data  from  an  insect  dispersal 
experiment,  in  which  the  movement  of  cabbage  root  flies  in  the  presence  of  a 
cabbage  crop  was  studied.  It  is  proposed  that  such  a  parameter  estimation 
method  can  be  a  useful  analytical  tool  to  help  develop  appropriate  models  in 
population  biology. 


I.  Introduction 


Vhen  constructing  a  mathematical  model,  one  often  finds  it  desirable 
to  estimate  unknown  parameters  in  that  model  based  on  actual  observations  of 
the  particular  phenomenon  under  study.  In  the  transport  equations  used  by 
mathematical  ecologists  to  model  population  dispersal,  one  typically  seeks 
estimates  of  parameters  representing  diffusion,  convection  and  "growth/death" 
rates.  Recognizing  that  these  models  are  very  naturally  posed  on  a 
two-dimensional  domain,  we  develop  a  bicubic  spline  based  parameter  estimation 
technique  for  a  general  two-dimensional  transport  equation.  Our  efforts  are 
guided  by  two  important  considerations:  one,  the  necessity  for  a 
computationally  efficient  procedure  and  two,  a  guarantee  that  such  a  procedure 
will  yield  a  set  of  parameter  estimates  that  is  "optimal"  in  some  sense. 

In  Section  II  we  give  a  precise  mathematical  formulation  of  what  we 
mean  by  an  estimation  problem.  It  will  be  seen  that  such  a  problem  is 
infinite  dimensional  in  nature  and  thus  difficult  to  solve;  so  one  is  led 
quite  naturally  to  define  a  sequence  of  finite  dimensional  "approximate 
estimation  problems"  in  an  attempt  to  obtain  a  solution.  A  theorem  is 
presented  relating  this  sequence  of  approximations  to  the  full  estimation 
problem,  and  the  result  serves  to  clarify  the  sense  of  "optimality"  alluded  to 
above.  The  actual  implementation  of  the  approximation  technique  formulated  in 
Section  II  will  be  described  in  Section  III  and  numerical  examples  will  be 
given.  Finally,  in  Section  IV,  we  present  the  results  of  an  analysis  of  real 
biological  data. 

Throughout  this  report  we  shall  adopt  the  following  conventions. 

Differentiation  of  a  function,  u,  with  respect  to  t,  x,  and  y  will  be  denoted 

by  u  ,  u  ,  and  u  ,  respectively,  with  obvious  extensions  to  higher  order 
i  x  y 

derivatives.  The  symbol  0  will  stand  for  the  Kronecker  product.  Ve  will  let 


2 


2  2 

L  =  L  (52)  be  the  usual  Hilbert  space  of  square  integrable  functions  on  52,  a 

2 

bounded  domain  of  R  ,  with  norm  ||  •  | |q  and  inner  product  <*,*>.  Likewise, 

03  00 

L  =  L  (52)  will  be  the  usual  Banach  space  of  essentially  bounded  functions 

with  norm  1 1  *  |  !„•  We  will  denote  by  H*(52),  i  =  0,1,2,...,  the  Sobolev 

2  th 

spaces  with  norm  ||  •  | | ^  of  "functions”  in  L  t52)  whose  i  "derivative"  is 
2  i 

also  in  L  (52).  Furthermore  Hq(S),  i  =  1,2,...,  will  be  the  Sobolev  space  of 
elements  of  H*(S2)  whose  "derivatives"  of  order  up  to  i-1  vanish  on  352,  the 
boundary  of  52.  Finally,  we  will  often  use  the  terms  parameter  and  coefficient 
interchangeably,  as  well  as  estimation  and  identification. 

II.  Mathematical  formulation 

We  consider  the  initial  value-boundary  value  problem  (IV-BVP)  on 
52  =  (0,1]  x  (0,1]: 

ut  +  7* (Vu)  =  7* (D  x  7u)  +  au  +  f,  tc(0,T], 

(2.1)  u(0,x,y)  =  uQ(y(x,y)), 

u(t,x,y)  =  0  on  3S2, 

where  f  =  f(0, t,x,y)  and  D,  V,  a,  and  0  are  assumed  to  be  functions  of  t,x, 
and  y,  with  D  =  (D^Dj)  with  V  =  (v^^).  The  dependent  variable,  u  = 
u(t,x,y),  represents  the  population  density  of  a  species,  whose  dispersal  over 
a  two-dimensional  domain  is  assumed  to  result  from  an  innate  diffusive 
mechanism,  D,  and  a  convective  or  "directed  transport"  mechanism,  V.  The 
parameter  a  represents  a  general  "source/sink"  or  "growth/death"  term.  (For 
the  definitive  account  of  diffusion  models  in  ecology  we  cite  the  excellent 
book  by  Okubo  (13]).  Without  loss  of  generality  we  have  assumed  homogeneous 
Dirichlet  boundary  conditions,  since  any  non-homogeneous  boundary  conditions, 
with  possibly  unknown  parameters,  can  be  included  in  the  parameters  0  and  y 
via  a  standard  transformation.  In  fact,  if  we  had  the  boundary  condition  u  = 


3 


F  on  92,  then  we  could  define  a  new  dependent  variable  to  be  v  =  u  -  g(x,y)’F 

where  g(x,y)  =  1  -  xy(l  -  x)(l  -  y).  Substituting  this  transformation  into 

(2.1)  to  get  an  equation  for  v,  we  see  that  we  must  assume  the  inhomogeneous 
1  2 

term  F  is  C  in  t,  and  C  in  x  and  y.  Notice  however,  that  g(x,y)  =  1  on  52 
so  that  v(t,x,y)  =  0  on  32.  Finally  we  note  that  solutions  to  (2.1)  will  be 
considered  in  the  sense  of  distributions  so  that  we  can  use  the  abstract 
framework  of  Lions  [12]  and  invoke  weak  variational  theory  to  discuss  the 
well-posedness  of  this  IV-BVP. 

Let  q  be  a  "vector"  of  parameters  (D,V,a, 3, y).  Dropping  the  Kronecker 
product  symbol®  for  ease  of  notation,  we  define  the  bilinear  form 

L(q) :  Hj(2)  x  hJ(2)  -►  R  by: 


(2.2) 


L(q )(♦,♦)  =  <DV$,  Vtj/>  -  -  <a^,^> 


and  then  rewrite  (2.1)  in  its  weak  form  in  the  state  space  L  (2): 


(2.3) 


<ut,4>>  +  L(q) (u, $)  =  <f,4>>,  for  all  in  H0(2), 


u(0)  =  Uq(y) , 


2  . 


where  <•,*>  is  the  standard  L  inner  product  on  [0,1]  x  [0,1].  Note  that  the 

boundary  conditions  have  been  incorporated  into  the  equation  of  state. 

We  assume  throughout  that  q  belongs  to  some  set  Q  of  admissible 

1  2  2 

parameter  functions  that  is  compact  in  the  space:  X  =  {H  ([0,T];L  (2)))  x 
0  3 

{H  ((0,T)  x  2))  .  Furthermore  we  make  the  following  assumptions: 


(Al)  0  is  a  bounded  subset  of  [L  ( ( 0 , T ]  x  8)}  ,  with  the  bound  m2; 

(A2)  There  exists  m^  >  0  such  that  for  any  q  e  Q,  the  components  of  D 

satisfy  D^  >  m^,  i  =  1,2  for  all  t,  x,  and  y; 


1 hA4\Af  -/V 


mmm 
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(A3)  The  inhomogeneous  term  f  is  C*  in  all  of  its  arguments,  the  initial 

2 

data  Uq  belongs  to  H  (Q),  and  the  map  y  -»  Uq(y)  is  continuous  from 
L2(C)  to  H2  (S); 

(A4)  For  any  q  in  Q,  Dt  and  c  Qp  a  bounded  subset  of  L^JO.T]  x  B)  with 
the  bound  m^. 

Theorem  2.1s  Under  the  assumptions  (Al),  (A2),  and  (A3),  the  initial  value 
problem  (2.3)  has  a  unique  weak  solution  u(t,*,*;q)  £  Hq(B).  Furthermore 
this  weak  solution  depends  continuously  on  the  initial  data  Uq. 

Proof:  We  appeal  directly  to  the  weak  variational  theory  of  Lions,  pp. 

100-111,  [12].  /// 

With  the  state  equation  (2.3)  in  hand  we  now  assume  that  for  each  time 
^  in  (0,T],  i  =  1,...,P  we  are  given  a  matrix  A(t..)  of  observations  (taken  at 
the  m*n  locations  (x^,y-^) , . . . ,  (xm,yn)) .  Associated  with  each  A(t^)  is  a 
matrix  r(tjjq)  =  (u(tj ,x^ jy^jq))  of  model  based  "predictions",  obtained  by 
evaluating  the  q-dependent  solution  u  at  the  points  (x^.y^)  1  <  j  <  m  and  1  < 
k  <  n.  Using  this  information  we  seek  a  vector  of  parameters,  q,  that  will 
yield  the  "best  fit"  of  the  model  to  the  data.  Thus,  we  formulate  the 
parameter  estimation  problem  as  follows: 

P  9 

(2.4)  minimize  J(q)  =  5EH  II  ~  r(t.;q)  || 

i*l  ii 

over  all  q  belonging  to  the  set  of  admissible  parameter  functions,  0.  Here 
the  norm  is  taken  to  be  the  Euclidean  norm  in  the  space  Rmn  and  it  is  clear 
that  J(q)  is  a  pointwise  least  squares  cost  functional, 
j  An  unavoidable  difficulty  now  arises,  for  the  minimization  problem  we 

|  have  formulated  is  inherently  an  infinite  dimensional  one  and  thus  hard  to 

i 

i 


solve.  Indeed,  ve  must  always  contend  with  the  infinite  dimensionality  of  the 

state,  u(t,*,';q),  when  we  try  to  solve  (2.3);  and  if  we  also  consider 

variable  coefficients,  we  must  then  deal  with  a  minimization  over  an  infinite 

dimensional  set  of  parameters.  Because  we  usually  cannot  solve  equation  (2.3) 

exactly,  we  would  like  instead  to  compute  an  approximate  solution  to  (2.3)  and 

thus  generate  a  matrix  of  "approximate  predictions"  to  use  in  the  cost 

functional  J.  We  would  also  like  to  employ  a  similar  strategy  to  reduce  the 

infinite  dimensional  minimization  problem  to  a  sequence  of  finite  dimensional 

minimizations.  Since  approximation  schemes  for  these  two  cases  will  be 

essentially  independent  of  one  another  [2],  it  will  greatly  simplify  our 

notation  if  we  focus  first  on  the  approximation  of  the  state  space  and 

postpone  discussion  of  the  parameter  set  approximation. 

Galerkin  schemes  based  on  cubic  splines  are  already  well  developed  for 

approximating  solutions  to  transport  equations  on  a  one-dimensional  domain 

( 1 ] ,  [2],  [3],  ]4],  [5]  and  these  are  readily  extended  to  two-dimensional 

domains.  Following  the  standard  Galerkin  technique,  we  define  a  sequence  of 

N  2 

finite  dimensional  approximating  state  subspaces  H  c  L  ,  N  =  1,2,...,  with 
N  2  N 

P  :  L  -»  H  the  canonical  orthogonal  projection.  Furthermore,  we  assume  that 
HN  c  Hq(Q)  and  that: 

(B)  For  all  z  e  C2(2)  n  hJ(Q),  | |PNz-z  |  |q  and  | |W(PNz-z) |  |Q 

are  less  than  or  equal  to  e<N) { | |zxx | |q  +  |  lzyy  I  Iq)  where 
c(N)  ■»  0  as  N  ■»  •. 

We  remark  that  by  using  arguments  from  Chapter  6  of  (14]  and  pp.  17-18  of  [6], 
it  can  be  seen  that  a  finite  element  approximation  scheme  based  on  bicubic 

CD 

splines  satisfies  these  conditions.  In  addition,  the  following  L  convergence 
rate  can  be  shown  to  hold:  )  |PNz-z  1 <  (c/N2)  { |  |zxx  1 1  „  +  ||zyy||j> 


(mm 
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N  N 

Restricting  the  bilinear  form  (2.2)  to  H  x  H  ,  we  then  define  the 
N  N 

Galerkin  approximation,  u  =  u  (t,*,*;q)  to  be  the  solution  of: 


(2.5) 


<u*,*>  +  L(q)(uN,\|«)  =  <f,^>,  for  all  \p  in  HN, 
u  (0)  =  P  Uq. 


We  note  that  we  can  view  this  equation  as  the  result  of  projecting  equation 

N 

(2.3)  into  the  finite  dimensional  subspace  H  ,  and  as  such  (2.5)  inherits  the 


well-posedness  of  its  parent  equation. 


Theorem  2.2:  Under  the  assumptions  (Al),  (A2),  and  (A3),  the  initial  value 

problem  (2.5)  has  a  unique  weak  solution  u  (t,’,*;q),  belonging  to  H  (Q)  and 

depending  continuously  on  the  initial  data.  In  fact,  given  the  finite 

N 

dimensionality  of  the  equation,  u  is  a  strong  solution. 

Proof:  As  before,  we  appeal  directly  to  the  theory  of  Lions,  pp.  100-111, 

[ 12  J .  /// 


Having  defined  a  sequence  of  equations  approximating  (2.3),  we  can 
define  a  sequence  of  approximate  estimation  problems  as  follows: 

p 

(2.6)  minimize  JN(q)  =  2ZZ  I!  A(^>  -  uN(t.,x  ,y  ,q)  |  |2 

i=l  1  1  J  K 

over  all  q  belonging  to  the  set  of  admissible  parameter  functions,  Q.  We  note 
however,  that  (2.6)  is  still  an  infinite  dimensional  problem,  and  this  fact 
leads  us  to  discuss  the  approximation  of  the  parameter  set  Q. 

In  a  manner  similar  to  that  for  the  state  approximation  scheme  we 

M 

define  a  sequence  Q  of  finite  dimensional  sets  that  approximate  0  in  a 
certain  sense.  Recall  that  we  have  defined  X  to  be  the  space: 

{H^( ( 0 , T ] ; L2(C) )} 2  x  (H^((0,T)  x  5)}"^.  Then,  following  the  ideas  developed  in 


7 


[2],  we  let  QM,  M  =  1,2,...,  be  subsets  of  X  defined  by  =  iM(Q),  where  i*1: 

Q  c  X  ■*  X.  In  addition  we  assume: 

M 

(Cl)  The  mappings  i  :  Q  -»  X  are  continuous; 

M 

(C2)  For  each  q  e  Q,  i  (q)  -*  q  as  M  -*  ®,  uniformly  in  Q. 

We  observe  here  that  in  Section  III  we  shall  discuss  an  approximation  scheme 
for  Q  based  on  linear  splines,  and  it  can  be  shown  that  such  a  scheme 
satisfies  the  above  assumptions.  In  fact  we  will  take  the  mappings  i  to  be 
the  M**1  order  interpolating  linear  spline  operator.  That  is  for  q  c  Q: 

iM(q)  =SH  q(ti,Xj,yk)S.(t)\.(x)vk(y)  i,j,k  =  0,1,.. .,M 

where  6.,  X.,  and  v,  are  the  standard  linear  splines  defined  on  equal 
1  J  K 

partitions  of  [0,T],  [0,1],  and  [0,1]  with  mesh  sizes  T/M,  1/M,  and  1/M, 

M 

respectively  (see  [14]).  With  these  definitions  Q  can  be  characterized  as: 

0M  -  {Q!  10, TJ  x  S  -  R5  I  ,n  -  3C  VjkVjV  vhere  °nijk  C  W1' 
Here  the  sum  is  of  course  taken  over  i,j,k  =  0,1,...,  M  and  are 

appropriate  compact  subsets  of  R.  We  also  note  that  we  have  taken  the  same 
degree  of  approximation,  M,  in  each  of  the  variables  t,x,  and  y,  but  only  for 
illustrative  purposes.  For  complete  details  of  this  construction  see  [2]. 

Given  these  approximating  sets  Q  ,  we  can  now  define  a  sequence  of 

finite  dimensional  approximate  estimation  problems  by  restricting  the 

N  M 

minimization  of  J  in  (2.6)  to  be  over  the  set  of  all  q  in  Q  .  The  following 

theorem  relates  the  solutions  of  these  problems  to  a  solution  of  the  full 

estimation  problem  (2.4). 

Theorem  2.3:  Let  Q  be  compact  in  X.  Then  for  each  M  and  N,  the  minimization 

-N  * 

problem  (2,6)  has  a  solution  q^.  Furthermore,  there  exists  an  element  q  of 


•vjrv; 


r't*. 
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-N  * 

0  and  a  subsequence,  denoted  again  by  qM,  that  converges  to  q  in  X,  and  this 
limit,  q*,  is  a  solution  to  the  full  parameter  estimation  problem  (2.3). 

We  will  not  present  the  entire  proof  of  this  theorem  here,  but  rather  will 
outline  just  the  essential  arguments.  For  complete  details  see  [17]. 

Let  uN(-,q)  be  the  Galerkin  solution  to  equation  (2.5).  It  can  be 
shown  that  for  each  te(0,T],  the  mapping  q  -»  u^(t,q)  from  the  set  of 

CO 

admissible  parameter  functions,  Q,  to  L  (Q)  is  continuous.  This  immediately 

N  M 

implies  the  continuity  of  the  approximate  cost  functionals  J  .  The  sets  Q 

are  compact  since  each  is  the  image  of  the  compact  set  Q  under  the  continuous 
M 

map  i  (see  assumption  (Cl)  earlier).  Hence  for  each  M  and  N,  there  exists  a 

N  N  M  H 

minimizer  qN  of  J  over  Q  .  Furthermore,  by  definition  of  the  sets  Q  ,  there 

"M  -M  H  ~M 

exists  a  sequence  (qN)  c  Q  such  that  qN  =  i  (qN).  Compactness  of  Q  then 

*  ft 

implies  the  existence  of  a  subsequence,  which  we  write  again  as  qN,  that 

converges  to  some  q*  in  Q.  It  can  then  be  shown  that  the  corresponding 

-M  * 

subsequence  q^  also  converges  to  q  . 

Finally,  to  argue  that  q*  solves  the  full  estimation  problem  (2.4) 

requires  proving  the  fundamental  result  that  "convergence  of  any  sequence  of 

parameters  q  to  some  q  implies  convergence  of  j”(q  )  to  J(q  )  as  N  and  K  •* 

N 

Since  the  cost  functionals  J  and  J  represent  pointwise  least  squares 

fit-to-data  criteria,  we  see  that  we  must  guarantee  pointwise  convergence  of 

our  state  approximations.  That  is,  under  an  appropriate  topology  on  Q,  we 

K  *  N  K 

must  show  that  convergence  of  q  to  q  implies  u  (tjX^y^;  q  )  converges  to 
u  ( t ,  Xj ,  y^*  q  )  as  N  and  K  -»  ®,  for  each  of  the  data  points  (Xj,y^)  1<  j  <  m 
and  1  <  k  <  n.  Notice  that  we  are  not  asking  for  global  pointwise  convergence 
of  the  state  approximations.  Rather,  it  will  suffice  to  show  local  pointwise 
convergence,  and  then  piece  together  at  most  a  finite  number  (in  fact  m-n)  of 
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such  results.  The  arguments  for  this  convergence  rely  on  a  general  technique 
from  finite  element  theory  that  estimates  a  local  L*  norm  in  terms  of  a  global 
norm.  We  then  use  the  weak  variational  framework  of  our  approximation 
scheme  to  derive  inequalities  from  which  we  can  obtain  appropriate  estimates 
to  establish  convergence  in  this  stronger  norm.  For  another  approach  to 
pointwise  convergence  in  two-dimensional  domains  see  [11] . 
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III.  Numerical  implementation  and  examples 

We  turn  now  to  a  description  of  the  implementation  on  the  computer  of 
the  approximation  framework  formulated  in  Section  II.  For  transport  equations 
on  a  one-dimensional  domain  cubic  spline  based  Galerkin  schemes  have  been  used 
successfully  and  their  computational  efficiency  well-documented  (1],  [2],  {31, 
[4],  [5J .  We  shall  see  shortly  that  the  key  feature  of  a  bicubic  spline  based 
approximation  technique  for  two-dimensional  transport  equations  is  that  the 
necessary  computations  reduce  in  a  natural  way  to  those  computations  arising 
in  the  one-dimensional  problem.  For  notational  simplicity  we  will  illustrate 
this  feature  for  the  case  of  constant  parameters,  then  describe  the 
straightforward  modifications  needed  when  variable  coefficients  are 
considered. 

Recall  that  we  have  defined  the  Galerkin  approximation  to  our  original 
IV-BVP  to  be  the  solution  of: 

<u>  +  <0^  ,  V\} »>  =  <Vu^,W+>  +  <otu^,^>  +  <z  f , 

(3.1)  M  M  M 

u (0)  .  P  uQ,  for  all  *  in  H  , 

where  we  have  rewritten  (2.5)  using  the  definition  of  the  bilinear  form 

L(q).  We  take  HN  to  be  the  span  of  (0.  .  .  n  (where  (3,.  *  B.(x)B.(y)  are 

ij  i,j=u  lj  l  3 

the  bicubic  splines  that  arise  as  pairwise  products  of  cubic  splines  B., 

N 

corresponding  to  the  partition  A  =  {i/N}^_Q  on  [0,1J  and  satisfying  B^(0)  = 

N 

B^(l)  =  0.  Thus  we  can  write  the  Galerkin  approximation  u  as: 

(3.2)  uN(t)  =  2C  v^.(t)gij(x,y)  i,j  =  0,1,2 . N 

M 

where  we  have  suppressed  the  dependence  of  (3^  on  N.  Notice  that  u 
automatically  satisfies  the  homogeneous  Dirichlet  boundary  conditions. 
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,N 


With  this  choice  of  the  subspaces  H  ,  equation  (3.1)  is  equivalent  to: 


<u 


N 


N 


",0.  .>  +  <D,u  ,  (0.  . )  >  ♦  <D„u  ,  (0.  . )  >  = 
t'  13  1  x’v  13  x  2  y  13  y 


(3.3) 


ij  y 


<v.uN,o, .)  >  +  <v0uN,<e, .)  >  +  <omN,e.  .>  +  <f,e.  .> 
1  ’ v  i3  x  2  i3  y  13  13 


N, 


<u"(O),0...>  =  <u  , pij> 


for  i , j  =  0,1, ... ,  N. 
N 


Substituting  the  above  expression  for  u  ,  equation  (3.2),  into  (3.3),  carrying 
out  the  necessary  differentiations,  and  using  the  linearity  of  the  inner 


product,  we  arrive  at  a  finite  dimensional  system  of  ordinary  differential 


N 


equations  for  the  Fourier  coefficients,  w^(t): 


(3.4) 


{ AN  (g)  AN]VN  =  -(D1(BN®  an]wn  +  D2[AN®BN]VN) 

+  v1[CN(g>  AN]WN  +  v2{AN®CN]WN 

+  «[an<siaV  +  fn, 

Ian®aV(0)  =  en, 

N 

where  (recall  B^  =  Bi  for  each  i) 

a”  =  |  B. (x)B. (x)dx,  B^  =  |  BC (x)Bj (x)dx, 

AJ  Jq  1  J  AJ  Jq  1  J 

Cij  ■  Ej.  .  J^B.(x)Bj(y)Uo(r)<lyd*, 

and  J~Bi(x)B.j(y,f(B>dydx, 


and  [M$N]  represents  a  Kronecker  product  of  matrices.  It  is  important  to 


N  N  N  2 

note  that  the  entries  of  the  matrices  A  ,  B  ,  and  C  are  just  the  pairwise  L 


inner  products  on  (0,11  of  the  one-dimensional  cubic  spline  basis  elements  and 
their  derivatives;  furthermore,  the  local  support  properties  of  cubic  splines 
guarantee  a  very  nice  banded  structure  for  these  matrices,  in  fact  they  are 
hepta-diagonal.  We  also  remark  that  we  have  factored  the  parameters  D^,  v^, 
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and  «  out  of  the  inner  products  in  equation  (3.3),  something  we  could  not  do 

if  the  coefficients  were  spatially  dependent.  However,  it  will  be  seen  that 

the  computational  feature  ve  are  trying  to  illustrate  is  independent  of 

whether  or  not  the  coefficients  are  constant. 

At  this  point  we  could  revrite  (3.4)  in  matrix-vector  form,  but  the 

computational  cost  of  numerically  integrating  the  resulting  system  would  be 

2 

prohibitive,  due  both  to  the  size  (dimension  =  (N+l)  )  and  the  sparse  nature 
of  the  matrices  involved.  Fortunately  however,  the  action  of  the  Kronecker 
product  [MQN]  on  any  matrix  Z  turns  out  to  be  given  by: 

(3.5)  lM<g)N]Z  =  M-Z-N* 


where  the  multiplication  on  the  right  hand  side  is  just  ordinary  matrix 
multiplication.  Hence  the  numerical  solution  of  (3.4)  requires  operation 
primarily  involving  the  matrices  that  arise  in  a  one-dimensional  cubic  spline 
approximation  scheme,  namely  AN,  B^,  and  C^.  Indeed,  we  have: 

an-wn.an  .  -(d/.«n.am  ♦  d/.^.b") 

(3'6)  +  VjC^.a"  ♦  v/.vNcV  ♦  «AN-VN.AN  ♦  FN, 

an.vn(0).an  .  en, 


bj  llj 

where  we  have  used  the  fact  that  A  and  B  are  symmetric.  It  is  here  that  the 
computational  attractiveness  of  this  scheme  really  becomes  apparent  because 
for  a  given  level  of  approximation  N,  ve  compute  and  store  ahead  of  time  the 
matrices  A^,  B^,  C^,  E^,  and  F^  in  eithei  band  or  band  symmetric  mode.  Then 
we  call  each  one  up  as  necessary  when  ve  are  solving  (3.6),  taking  advantage 
of  their  banded  structure  to  speed  up  the  matrix  multiplications.  Once  we 
solve  equation  (3.6)  (equivalently  (3.4)),  we  can  reconstruct  the  approximate 
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solution  (3.2)  and  generate  a  set  of  "predicted  observations"  to  be  compared 
with  the  real  data.  However,  before  ve  discuss  the  computer  code  that 
actually  integrates  (3.6)  and  produces  the  Galerkin  approximation  used  in  the 
minimization  of  the  cost  functional  J^(q),  we  describe  the  modifications  to 
(3.4)  and  (3.6)  that  must  be  made  in  the  presence  of  variable  coefficients 
12). 

If  we  include  only  temporal  variation  in  the  parameters,  there  is 
little  to  change,  since  we  can  still  factor  the  coefficients  out  of  the  inner 
products  in  equation  (3.3).  Taking  a  linear  spline  representation  for  each 
of  the  parameters,  e.g.  D^(t)  = ^  6jL(t),  k  =  1,2,...,M,  we  simply 
replace  Dp  Dp  Vp  v«,  a,  0,  and  y  by  their  appropriate  representations,  thus 
making  (3.6)  a  non-autonomous  equation  due  to  the  linear  spline  functions 
lk(t).  Having  solved  this  time-dependent  analogue  to  equation  (3.6)  and 
reconstructed  the  approximate  solution  (3.2),  we  then  seek  to  minimize  the 
cost  functional  JN(*)  over  some  finite  dimensional  set  of  Euclidean  parameters 
containing  the  admissible  values  for  the  coefficients  in  the  various  linear 
spline  representations  for  Dp  Dp  Vp  Vp  o,  0,  and  y. 

Turning  to  the  case  of  spatially  varying  parameters  we  see  that  it  is 

not  quite  so  simple  to  handle  because  we  can  now  no  longer  factor  the 

coefficients  directly  out  of  the  inner  products  in  (3.3).  Nevertheless  we  can 

preserve  the  Kronecker  product  structure  and  action  given  by  (3.5).  More 

importantly,  we  will  still  be  able  to  compute  and  store  ahead  of  time  a  set  of 

N  N  N 

"inner  product  matrices",  analogous  to  A  ,  B  ,  and  C  from  (3.4)  to  be  used  in 
the  solution  of  the  counterpart  to  the  system  (3.6). 

As  before  we  take  a  linear  spline  representation  for  each  of  the 
parameters,  this  time  in  x  and  y.  To  illustrate,  we  examine  what  happens 
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N 


to  the  first  term  <D,u  , (0..)  >  of  equation  (3.3).  For  ease  of  notation 

IX  ij  x 


we  assume  our  coefficients  are  separable;  however  we  hasten  to  aid  that  the 
method  will  still  apply  in  the  more  general  case  of  non-separable  parameter 
functions.  Ve  write: 


D^x.y)  =  D^(x) *D^(y)  =  (Eyk<*>>  '  (Eyk(y>>  k  =  1,2,...,M. 

N 


Substituting  in  this  expression  for  and  the  representation  (3.2)  for  u 
differentiating  where  appropriate,  invoking  the  linearity  of  the  inner 
product,  and  using  the  separability  of  the  coefficient  D^(x,y),  we  have  that: 

{<DlU”,(0i;.)x>)  =  [BN(DX)<£  AN(Dy)]VN 
=  BN(D^)-WN-(AN(D5[))t, 

where  A^D*)  =  E\  J^Bi(z)Bj(z)lk(z)d2  k  =  1,2 . M 

and  B?.(DX)  =  £>k  J^B'(z)Bj(z)lk(z)dz  k  =  1,2 . M. 

In  a  similar  manner  we  can  analyze  the  changes  necessary  in  the  remaining 
homogeneous  terms  of  equation  (3.3)  when  the  appropriate  linear  spline 
representation  is  substituted  for  the  parameter  appearing  in  each.  In 
particular  we  have  that: 


{<D2Uy,(0ij)y>}  =  IAN(D2)  ©bV^JV*1 


=  AN(D£)  •  V*  •  (BN(D2))t , 


{<vlUN,(0  )x»  =  lCN(v^)  ®  AN(v][)IVN 


=  CN(vX)  •  /  •  (AN(v5[))t, 


(<v2uN,(0ij)y>)  =  (AN(vX)©  CN(Vy)JWN 

=  AN(v2 )  •  WN  •  (CN(v^) ) * , 


and 


{<auN,0ij»  =  IAN(oX)©AN(ay))VN 


=  AN(aX)  •  WN  •  (A N(«y))t. 


N  N  N 

Here  the  matrices  A  (•)>  B  (•)>  and  C  (•)  are  just  linear  combinations  of  the 


elementary  matrices  A^,  B^,  and  C^,  whose  entries  are,  in  analogy  to  those  of 

2 

(3.4),  the  pairwise  L  inner  products  on  [0,1]  of  the  one-dimensional  cubic 
spline  basis  elements  and  their  derivatives,  now  weighted  by  the  k1*1  linear 
spline,  lk(*)»  k  =  1,2,. ..,H. 

We  note  that  it  is  the  convective  terms  of  (3.3)  that  give  rise  to  the 
matrices: 

cV)  =  {E\J^B.(z)Bj(2)lk(z)dz}  k  =  1,2 . M. 

In  addition,  we  observe  that  the  local  support  properties  of  cubic  splines 

again  guarantee  a  nice  hepta-diagonal  structure  to  these  "weighted  inner 

product  matrices";  furthermore  for  a  given  level  of  approximation,  N  for  the 

state  space  and  M  for  the  parameter  space,  we  can  compute  and  store  these 

matrices  ahead  of  time  just  as  in  the  simpler  constant  coefficient  case 

described  earlier.  We  point  out,  however,  that  if  the  inhomogeneous  term  f(0) 

and  the  initial  condition  u  (y)  are  non-linear  functions  of  the  parameters  £ 

N 

and  r,  then  their  respective  projections  into  the  subspace  H  (i.e.  the 
N  N 

matrices  E  and  F  of  (3.6))  cannot  be  computed  and  stored  ahead  of  time,  but 
rather  must  be  continually  updated  during  the  minimization  procedure. 
Nevertheless  this  needs  to  be  done  only  at  the  start  of  each  integration  of 
the  approximating  system  of  ordinary  differential  equations,  and  the  resulting 
matrices  will  still  exhibit  a  banded  structure.  Finally,  we  note  that  if 
temporally  and  spatially  dependent  parameters  are  considered  simultaneously, 
our  method  still  works  since  we  would  take  a  general  representation: 

q(t,x,y)  *  <Ebkyx>>  '  (  Eckyy))  k  = 

and  then  immediately  factor  the  temporal  dependence  of  q  out  of  the  inner 
product  and  proceed  as  before. 


We  now  discuss  the  actual  solution  of  the  ordinary  differential 

equation  (3.6).  Since  it  is  well  known  that  approximation  schemes  for 

parabolic  partial  differential  equations  often  give  rise  to  systems  of  stiff 

ordinary  differential  equations,  we  use  the  IMSL  code,  DGEAR,  with  its  option 

for  stiff  systems.  As  a  part  of  this  solution  process  the  general  equation 

[M®N]Z  =  Y,  equivalently  M-Z^N)1  =  Y,  must  be  solved  for  Z.  This  task  is 

accomplished  by  two  successive  applications  of  a  Cholesky  decomposition 

routine  that  enables  one  to  solve  FX  =  G  for  X.  In  fact  one  first  lets  X  = 

Z*(N)*  and  solves  MX  =  Y  by  a  Cholesky  decomposition  of  M  followed  by 

back-substitution.  Then  one  solves  N-fZ)11  =  X1  by  a  Cholesky  decomposition  of 

N  again  followed  by  back-substitution.  The  minimization  of  the  cost 
N 

functional  J  (q)  is  performed  by  another  standard  IMSL  code,  called  ZXSSQ, 
which  employs  a  modified  Levenberg-Marquardt  algorithm.  The  FORTRAN  code 
built  around  these  two  large  IMSL  routines  is  a  modification  of  an  original 
program  written  by  James  Crowley  [8].  The  testing  of  this  code,  which  we 
shall  now  describe,  was  performed  on  the  computing  system  at  Brown  University 
running  either  an  IBM  370  or  IBM  3081. 

In  the  following  examples,  we  considered  a  number  of  "data"  sets 
generated  from  various  known  solutions  of  the  constant  coefficient  model  on 
Q  =  [0,11  x  [0,1]: 

ut  -  D(uxx  *  V  *  Vx  *  v2uy  *  t>0 

(3.7)  u(0,x,y)  =  uQ(x,y) 

u(t,x,y)  =0  on  32 

Values  for  D,  v^,  and  «  were  specified,  then  a  simple  separation  of 
variables  technique  was  used  to  calculate  each  explicit  solution.  At  each  of 
the  times  t  =  .2,  .4,  and  .6  a  matrix  of  "observations"  was  then  obtained  by 


evaluating  these  solutions  at  the  (x,y)  grid  points  x  =  .1,  .3,  .5,  .7,  .9  and 
y  *  .1,  .3,  .5,  .7,  .9. 

We  carried  out  a  series  of  tests  on  each  data  set  in  which  various 

combinations  of  parameters  were  considered  to  be  unknown  and  the  remaining 

parameters  held  fixed  at  their  true  values.  Initial  guesses  of  the  "unknown" 

parameters  were  input  into  the  estimation  code  and  a  particular  level  of 

approximation  specified.  The  performance  of  the  code  could  then  be  assessed 

in  light  of  factors  such  as  the  "error"  made  in  the  initial  guess,  the  degree 

of  approximation  chosen,  or  how  many  parameters  were  considered  unknown.  We 

note  that  low  levels  for  the  state  approximations  were  taken,  N  =  4,6,  or  8. 

This  was  done  because  the  computational  cost  for  any  given  two-dimensional 

estimation  problem  far  exceeds  the  cost  for  a  comparable  one-dimensional 

problem  (see  results  of  numerical  experiments  in  [1],  [2],  [3],  and  [5]).  Of 

course  the  reason  for  this  is  obvious;  namely  that  the  two-dimensional  setting 

2 

gives  rise  to  a  system  (3.6),  that  consists  effectively  of  (N  +  1)  ordinary 
differential  equations,  as  opposed  to  (N  +  1)  in  the  one-dimensional  case.  In 
fact,  just  to  integrate  equation  (3.6)  for  N  =  12  takes  roughly  12  minutes  of 
CPU  time  on  the  IBM  370  and  for  N  =  14  that  figure  increases  to  40  minutes. 

If  the  faster  IBM  3081  is  used,  these  times  decrease  to  3  and  10  minutes,  but 
are  still  very  large  compared  to  a  one-dimensional  problem.  Thus,  with  an  eye 
towards  future  real-time  applications  of  our  estimation  technique,  we  believe 
it  would  be  fruitful  to  pursue  the  implementation  of  our  code  on  vector 
machines  and  other  supercomputing  processors. 

We  also  remark  that  the  efficiency  of  the  Levenberg-Marquardt 
minimization  algorithm  decreases  as  the  number  of  unknown  parameters 
increases.  This  is  particularly  significant  if  we  want  to  estimate  variable 
coefficients.  Indeed,  a  simple  constant  coefficient  model  such  as  (3.7)  has 
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oi^ly  four  unknown  parameters.  However,  if  we  take  a  general  transport 

equation  such  as  (2.1),  and  implement  the  parameter  approximation  scheme  we 

have  described,  then  the  dimension  of  the  unknown,  and  in  this  case 

approximate,  parameter  space  will  increase  dramatically.  For  example,  if  we 

consider  a  constant  diffusion  model  with  temporally  dependent  convection  and 

"growth/death"  terms  V^(t),  Vjft),  and  a(t),  and  we  assume  a  representation 

for  each  as  a  linear  combination  of  four  linear  splines,  then  we  will  have 

thirteen  unknown  parameters  to  search  for.  Obviously  this  figure  will 

increase  if  we  take  a  linear  combination  of  more  than  four  splines  in  our 

representations;  and  if  we  also  include  spatial  dependence  in  the  convection 

and  "growth/death"  terms,  then  we  must  deal  with  an  even  larger  number  of 

N 

unknown  parameters.  But  the  difficulty  of  minimizing  J  (•)  over  a  large  set 
of  unknowns  is  not  insurmountable,  especially  in  light  of  recent  advances  in 
computer  technology.  In  fact,  as  new  software  (different  minimization 
routines)  and  hardware  (array  and  parallel  processors)  become  available,  it 
is  likely  that  this  aspect  of  this  problem  can  be  treated  in  a  relatively 
efficient  manner. 


Example  3.1  We  considered  the  standard  heat  equation  with  the  diffusion 
coefficient,  D,  equal  to  1.0  and  the  initial  condition  given  by  the  product  of 
two  one-dimensional  "hat"  functions.  Convergence  results  were  quite  good,  see 
Table  3.1,  with  no  significant  improvement  in  accuracy  obtained  by  increasing 
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Example  3.2  In  this  example  we  chose  D  =  1.0,  v^  =  -2.0,  Vj  =  -2.0,  and  a  = 
19.7392  in  equation  (3.7)  with  uQ(x,y)  =  exsinnxe-’rsinity.  In  the  tests 
summarized  in  Table  3.2  the  level  of  the  state  approximation  was  held  at  N  =  4 
and  it  is  seen  again  that  the  method  produced  good  results. 

Example  3.3  For  this  example  we  set  D  =  .20,  v^  =  -1.0,  V£  =  -1.8,  and  a  = 
7.2478.  The  relative  magnitudes  of  the  coefficients  were  chosen  to  reflect  a 
typical  example  of  population  dispersal;  and  to  ascertain  whether  our 
technique  would  be  able  to  estimate  an  asymmetric  model  (which  is  almost 
always  the  case  with  real  biological  data),  a  directionally  dependent 
convective  component,  v^  not  equal  to  V£,  was  considered.  In  Table  3.3  we 
present  a  comparison  of  the  method's  performance  at  two  levels  of  state 
approximation  and  we  see  that  the  parameter  estimates  were  almost  all  within 
one-tenth  of  one  percent  of  the  true  values.  Thus,  an  improved  fit-to-data  is 
obtained  at  the  expense  of  increased  CPU  time,  but  with  no  great  improvement 
in  parameter  estimates. 

Example  3.4  In  this  example  we  considered  the  problem  of  estimating 
parameters  in  the  presence  of  "noisy"  data.  We  took,  the  same  model  as  in 
Example  3.3  but  introduced  random  error  into  the  "observations"  (see  Appendix 
for  details).  Representative  results  are  displayed  in  Table  3.4  and  Table 
3.5,  and  together  with  other  tests  these  suggest  that  our  technique  can 
perform  well  given  data  containing  noise.  It  is  an  interesting  qualitative 
result  that  in  test  3  of  Table  3.5  the  method  m?  >ged  to  estimate  the  correct 
relative  magnitudes  of  the  true  parameters. 
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Table  3.4 

Initial 

guess 

Estimated 

value 

True 

value 

RSSQ 

Iterations 

CPU(IBM  3081) 

1) 

D° 

ss 

.10 

D4  = 

.2003 

D  =  .20 

.6835 

5 

45  sec 

2) 

D° 

= 

.10 

D4- 

.20057 

D  =  .20 

.479 

9 

90  sec 

= 

-.50 

<i 

I-*  4^ 

II 

-1.0031 

v  x  =  -1.0 

v2 
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-4 

V2  = 
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N  a  4 


24 


Table  3.5 
(Noisy  data) 
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II 

13.4234 

o  =  7.2478 

N  =  4 


IV.  Parameter  estimation  in  a  two-dimensional  model  of  insect  dispersal 

In  a  series  of  mark-recapture  experiments,  Kareiva  recorded  the 
movement  of  flea  beetles  along  linear  arrays  of  collard  patches  [10].  A  cubic 
spline  based  parameter  identification  technique  was  subsequently  used  to 
analyze  this  data  and  to  assess  the  appropriateness  of  various  forms  of  a 
general  transport  equation  as  models  for  the  observed  movement  [3],  [4],  [5]. 
While  this  analysis  proved  quite  successful  in  identifying  different 
mechanisms  of  dispersal  and  quantifying  their  relative  importance,  it  is 
important  to  note  that  the  experimental  design  restricted  movement  to  a  one 
dimensional  domain,  thus  allowing  consideration  of  only  one-dimensional 
transport  equations  as  models.  It  is  apparent  of  course  that  a 
two-dimensional  domain  provides  a  more  natural  setting  for  most  models  of 
insect  dispersal  (indeed  for  most  models  of  population  dispersal)  and  in  this 
section  we  describe  the  application  of  our  estimation  technique  to  the 
analysis  of  cabbage  root  fly  dispersal  on  a  two-dimensional  domain. 

Our  data  are  taken  from  mark-recapture  experiments  by  Hawkes  in  which 
cabbage  root  flies  (Erioischia  brassicae)  were  related  at  a  point  adjacent  to 
and  downwind  from  a  cabbage  (brassica)  crop  [9].  Although  Wright  [16]  had 
rejected  anemotaxis  as  a  mechanism  of  attraction  and  Thornsteinson  [15] 
claimed  "there  seems  to  be  no  critical  evidence  that  insects  orient  to  plants 
beyond  a  few  meters",  wind  tunnel  experiments  by  Coaker  and  Smith  [7] 
indicated  that  female  E.  brassicae  do  fly  upwind  in  the  presence  of  brassica 
odor.  To  resolve  this  issue  Hawkes  sought  to  calculate  dispersal  rates  of  E. 
Brassicae  released  from  a  point  exposed  to  brassica  odor. 

When  recapture  data  suggested  random  dispersal,  an  empirical  model, 

In  y  =  a-  bTx,  was  used  to  relate  the  number  of  flies  captured,  y,  to  the 
distance  from  the  release  point,  x,  where  a  and  b  are  constants.  Average 
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distances  dispersed  and  thus  average  rates  of  dispersal  were  then  calculated 
on  the  basis  of  this  model.  However,  in  the  case  of  non-random  movement,  no 
such  model  was  available  to  estimate  rates  of  dispersal.  Ve  have  subsequently 
applied  our  parameter  estimation  technique  to  estimate  dispersal  rates  for 
this  type  of  recapture  data,  focusing  on  Hawkes'  data  for  gravid  (i.e.  egg 
bearing)  female  £.  brassicae,  as  this  group  was  observed  to  exhibit  the 
greatest  non-random  movement.  (We  are  quick  to  point  out,  however,  that  the 
applicability  of  this  estimation  technique  does  not  depend  on  any  a  priori 
knowledge  of  a  population's  specific  behavior.  Indeed, it  is  precisely  this 
behavior  that  we  usually  wish  to  ascertain  and  then  analyze). 

The  experiments  were  carried  out  in  a  large  field  bordered  on  the 
north  by  a  hedge.  A  30  by  30  meter  cabbage  plot  was  planted  immediately  south 
of  the  hedge  with  large  areas  of  fallow  ground  to  the  south,  east,  and  west  of 
the  plot.  Water  traps  spaced  six  meters  apart  were  placed  along  the  hedge, 
within  the  crop,  and  in  the  surrounding  fallow  area  as  shown  in  Figure  4.1. 
Since  the  prevailing  wind  direction  was  from  the  east-southeast,  flies  were 
released  from  a  point  at  the  hedge  24  meters  to  the  west  of  the  northwest 
corner  of  the  cabbage  crop.  Direct  observation  and  recapture  data  showed  that 
movement  did  not  begin  until  29.5  hours  following  the  time  of  initial  release. 
After  the  onset  of  dispersal,  data  representing  the  distribution  of  the  flies 
was  collected  during  two  consecutive  seven  hour  periods.  (For  further  details 
of  the  experiment  see  Hawkes'  paper  l 9 ] ) . 

We  considered  the  following  two  dimensional  transport  equation  as  a 
model  for  describing  the  distribution  of  the  flies: 

(4.1)  ut  =  D(uxx  +  uyy)  -  vx(t)ux  -  v2(t)uy  -  oc(t)u 

(x,y)  c  (0,1]  x  (0,1], 


t  >  0. 
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Initial  data  were  given  by  u(0,x,y)  =  0  for  (x,y)  *  (.5, .5)  and  u(0,.5,.5)  = 
1930  (the  number  of  marked  flies  released).  Since  there  were  no  cabbage 
plants  outside  the  furthest  extent  of  Hawkes'  trap  grid,  we  assumed  Dirichlet 
boundary  conditions,  u(t,0,y)  =  u(t,l,y)  =  u(t,x,0)  =  u(t,x,l)  *  0.  Here  the 
variables  x  and  y  represent  dimensionless  quantities  based  on  the  scaling: 
x  =  x  meters/127.5  meters  and  y  =  y  meters/127.5  meters.  Thus  the  entire 
field  was  rescaled  to  fit  on  the  unit  square,  with  the  release  point 
corresponding  to  x  =  y  =  .5.  The  westernmost  traps  and  the  easternmost  traps 
corresponded  to  y  =  .3353  and  y  =  .9,  respectively,  while  the  northernmost 
traps  (those  along  the  hedge)  and  the  southernmost  traps  corresponded  to  x  = 

.5  and  x  =  .6647.  Time  was  rescaled  as  t  =  t  hours/24  hours,  so  that  tj  = 
.14583  corresponded  to  the  midpoint  of  the  first  seven  hour  census  period 
following  the  beginning  of  dispersal  and  t ^  *  .4375  corresponded  to  the 
midpoint  of  the  second  seven  hour  census  period.  The  convection  coefficients, 
Vj(t)  and  V2(t),  as  well  as  the  "growth/death"  term,  «(t),  were  represented  by 
a  linear  combination  of  four  linear  splines.  The  function  o(  t )  was  assumed 
positive  so  that  the  entire  negative  term  a(t)u  would  reflect  the 
disappearance  of  flies  from  the  experiment  (e.g.  through  actual  death,  long 
range  migration,  wearing  off  of  the  radioactive  marker,  etc.). 

Ve  began  our  analysis  by  considering  the  simpler,  constant  coefficient 
equation  as  a  model.  Then,  by  allowing  one  or  more  of  the  coefficients,  v^, 
V2»  and  a,  to  vary  in  time,  we  gradually  increased  the  complexity  of  the  model 
to  arrive  at  the  general  equation  (4.1).  At  each  stage  we  sought  to  minimize 
the  residual  sum  of  squares  of  differences  (RSSO)  between  model  predictions 
and  observed  data.  Although  it  is  true  that  ve  could  have  considered  the  full 
equation  (4.1)  from  the  onset  and  asked  the  computer  to  identify  all  the 
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parameters  in  the  model  at  once,  ve  emphasize  that  this  is  not  the  procedure 

ve  followed  (see  our  earlier  remarks  in  Section  III  regarding  the  minimization 

algorithm  employed).  Instead,  ve  estimated  various  combinations  of  parameters 

while  holding  the  remaining  coefficients  fixed  at  certain  "nominal"  values. 

As  a  rule,  these  were  taken  to  be  the  "best"  values  returned  from  the 

consideration  of  a  simpler  model.  For  example,  to  identify  the  time-dependent 

profile  of  the  "death"  term  o(t),  ve  first  fixed  the  diffusion  and  convection 

parameters  at  those  values  that  provided  the  "best  fit"  of  the  constant 

coefficient  model  to  the  recapture  data,  and  then  we  estimated  the  "weighting 

3 

factors",  cl,  in  the  linear  spline  representation  a(t)  =  ^T]a._l,_(t).  All  of 

k=0 

the  computational  work  for  this  analysis  was  performed  on  the  CDC  6600 
computer  at  Southern  Methodist  University  and  the  IBM  3081  at  Brown 
University. 

In  the  constant  parameter  model  we  sought  to  estimate  the  four 
parameters  D,  v^,  v^,  and  a.  This  analysis  immediately  revealed  a  number  of 
interesting  qualitative  features  (see  Table  4.1): 

(1)  a  relatively  small  diffusive  term,  D; 

(2)  the  presence  of  some  directional  bias  in  the  convection  terms, 

vi  <  v 

(3)  the  importance  of  a  large  "death"  term,  a. 

These  results  are  not  surprising,  in  fact  they  correlate  quite  well  with  what 
even  a  casual  perusal  of  the  data  suggests  (see  Figure  4.1).  There  is  little 
evidence  of  purely  random  dispersal,  indeed  some  directed  movement  is  almost 
surely  present.  Also  the  fact  that  only  fifty-five  out  of  the  original  1930 
flies  released  were  recaptured  at  the  first  census  certainly  indicates  the 
need  for  a  large  "death"  term. 


We  next  introduced  temporal  variation  in  the  coefficients  beginning 

3 

with  the  "death"  term,  a(t)  =  ]JT]  «.l.(t),  where  l.(t)  are  the  standard 

k=0 

linear  splines  as  in  Sections  II  and  III.  The  terms  D,  v^,  and  v2  were  held 

fixed  and  we  sought  estimates  of  the  coefficients,  o^,  that  would  lead  to  a 

reduction  in  the  RSSQ  from  the  constant  parameter  model  (see  Table  4.2). 

Though  we  were  able  to  obtain  some  profiles  for  a(t)  that  could  be  considered 

biologically  plausible,  we  could  not  reduce  the  RSSQ  significantly  (RSSQ  =  200 

for  the  constant  coefficient  model,  RSSQ  =  198  for  the  variable  a  model). 

3 

Time  varying  convection  terms  were  considered  next,  v.(t)  =  &lr.(t) 

3  1  k=0  K  K 

and  v„(t)  =2]  YjIk ( t ) .  But  as  before,  though  we  obtained  some  bio- 
k=0  K  K 

logically  arguable  profiles  for  the  convection  terms  (see  Table  4.3),  we  could 
not  reduce  the  RSSQ.  What  we  did  observe  consistently  was  a  general 
intensitivity  of  the  RSSQ  to  the  latter  two  coefficients  in  each 
representation  for  v-^t),  v2(t),  and  ot(t).  This  indicated  that  our  models 
were  simply  not  predicting  the  data  at  the  second  time  point.  At  this  stage 
it  was  thought  that  the  diffusion  coefficient , though  already  very  small,  might 
nonetheless  also  exhibit  some  temporal  dependence.  Accordingly,  some  tests 
were  performed  to  estimate  a  time  varying  diffusion  term,  but  still  there  was 
no  corresponding  improvement  in  the  model's  prediction  of  the  data  (RSSQ 
remained  at  198).  In  fact,  it  is  interesting  to  note  that  for  this  data  set, 
the  total  sum  of  squares,  TSSQ,  equals  161.  This  implies  that  simply  using 
normally  distributed  noise  as  a  model  produced  a  better  fit-to-data  than  any 
of  the  "dynamic"  models  we  had  considered  so  far. 

However,  these  tests  were  not  without  value,  for  they  convinced  us 
that  our  difficulties  lay  with  the  magnitude  and  shape  of  the  variable  "death" 
term,  a(t).  As  a  remedy  we  introduced  a  discontinuous  a(t)  into  the  model 
with  the  discontinuity  located  at  the  first  time  point  t^.  We  retained  the 


profile  of  a(t)  (taken  from  earlier  testing)  from  t^  to  t^,  then  set  a(t)  s  0 
for  t  >  tjj  in  effect  "turning  off"  the  "death"  mechanism  after  t^.  The 
biological  interpretation  of  such  a  profile  is  that  immediately  following  the 
time  of  release  there  is  a  large  emigration  of  the  population  due  to  the 
disturbance  of  marking  and  handling  ( 5 ] .  Emigration  then  settles  down  to  its 
more  "intrinsic"  low  level,  which  we  artifically  took  to  be  zero.  It  was 
hoped  that  this  would  allow  the  convective  and  diffusive  mechanisms  to 
redistribute  the  population  over  the  next  time  period.  In  fact,  the 
introduction  of  this  discontinuity  resulted  in  a  significant  decrease  n  the 
RSSQ. 

Since  we  had  already  determined  the  purely  diffusive  component  of 
dispersal  to  be  small,  we  held  the  diffusion  term  fixed  at  its  "best" 
estimated  value  from  the  constant  coefficient  model  testing.  Ve  also  held  the 
"death"  function  fixed  at  the  discontinuous  profile  described  above.  By 
estimating  various  combinations  of  the  coefficients  3^  and  y^  in  the 
representations  of  the  terms  v^  and  V£,  we  were  able  to  produce  convection 
profiles  that  reduced  the  RSSQ  to  a  value  of  98  (see  Table  4.4).  In  light  of 
this  result,  the  insensitivity  of  the  RSSQ  to  o^,  and  y^,  i  =  2,3,  in  the 
expressions  for  v^,  V2,  and  a  can  be  understood.  Before  the  introduction  of 
the  discontinuous  "death"  term,  the  very  large  value  of  a(t)  at  t^  influenced 
the  model's  behavior  for  a  significant  amount  of  time  after  t^  (even  if  oc,  * 

=  0.0).  This  large  value  dominated  the  mechanisms  of  dispersal  and  caused 
the  model  to  predict  a  population  identically  equal  to  zero  at  the  second  time 
point,  tj,  regardless  of  the  profiles  of  the  convection  terms.  Only  by 
"turning  off  the  decay",  via  the  discontinuity,  could  we  allow  for  convection 
to  be  identified  as  a  significant  component  of  the  motion  (see  Figure  4.2). 
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Having  produced  a  model  which  allowed  more  freedom  to  the  dispersive 

mechanisms,  we  proceeded  to  carry  out  a  set  of  tests  to  try  to  "fine  tune"  the 

convective  terms.  Since  some  of  the  profiles  that  we  had  estimated  yielded  a 

good  RSSQ  but  made  little  biological  sense,  we  sought  to  estimate  more 

biologically  reasonable  convection  functions.  The  results  depicted  in  Table 

4.5  show  that  our  attempts  met  with  qualified  success.  While  the  qualitative 

feature  of  decreasing  x-convection  anu  increasing  y-convection  that  we 

identified  can  be  explained  by  the  greater  proximity  of  the  cabbage  crop  to 

the  release  point  when  measured  in  the  x-direction  than  when  measured  in  the 

y-direction,  we  are  quick  to  note  that  this  set  of  profiles  did  not 

significantly  reduce  the  RSSQ  from  its  previous  best  minimum  value  of  98.  In 

fact,  we  actually  identified  several  sets  of  convection  functions  during  our 

analysis  that  produced  this  value  for  the  RSSQ.  Such  an  example  of 

non-uniqueness  is  of  course  not  unexpected,  being  a  reflection  of  the  inherent 

"ill-posedness"  of  many  inverse  problems  of  this  type. 

We  also  performed  tests  to  try  to  "fine  tune"  the  variable  "death" 

term.  Recognizing  the  artificial  nature  of  the  discontinuity  introduced  in 

a(t) ,  we  sought  to  identify  a  profile  for  a(t)  with  a  steep  gradient  in  a 

neighborhood  of  the  first  time  point  t, .  Here  we  considered  the 

11 

representation  a(t)  =  o^l^ft)  and  set  otg  a  and  m  s  a  . . .  * 

k=0  1 

<x^.  We  then  estimated  the  "weighting  factors"  a^,  a^,  a^,  and  <x^.  The 
results  are  displayed  in  Table  4.6  and  we  see  the  presence  of  a  very  steep 
gradient  in  the  estimated  profile  for  a(t). 

Finally,  we  have  performed  a  series  of  tests  to  identify  spatial 
dependence  in  the  convection  terms.  We  assumed  the  representations 
3  3  3  3 

v^t.x)  -  (E  ekik<o)  (Eck1k(x))  and  v2(t>y)  =  (E  vkilt(t))-(E\1|C(y)) 
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where  the  coefficients  0^  and  Y^  were  set  identically  equal  to  1.0  so  that 
purely  spatial  variation  in  convection  was  considered.  Our  estimation 

and  see  Table 

4.7.  Next,  using  the  resulting  spatially  dependent  profiles,  new  estimates 

and  defining  the  temporally  dependent 
portions  of  convection,  see  Table  4.8.  We  caution,  however,  that  while  the 
inhomogeneity  of  the  experimental  site  certainly  suggests  that  we  consider 
spatial  dependence  in  the  convective  terms,  the  experimental  feature  of  a 
single  point  release  does  not  allow  the  proper  separation  of  temporal  effects 
on  convection  from  spatial  effects.  So  while  we  actually  succeeded  in 
lowering  the  RSSQ  to  92.42  (see  Table  4.7)  and  subsequently  to  89.68  (see 
Table  4.8),  it  is  not  clear  that  our  data  can  support  these  results. 

We  turn  now  to  a  discussion  of  a  statistical  criterion  which  we  used 
to  help  evaluate  the  relative  strengths  of  our  models.  Following  the  method 
described  in  [4]  and  [5]  (an  ad  hoc  modification  of  multiple  regression 
analyses  and  significance  tests  based  on  the  F-distribution) ,  we  calculated 
F-statistics  comparing  the  variation  explained  by  a  particular  model,  TSSQ  - 
RSSQ,  with  the  unexplained  variation,  RSSQ.  The  degrees  of  freedom  for  these 
two  quantities  were  taken  to  be  k  and  (n  -  k  -  1),  respectively,  where  k 
equals  the  number  of  unknown  parameters  in  the  given  model  and  n  equals  the 
number  of  data  points.  Hence,  recalling  that  the  TSSQ  for  our  data  set  is 
161,  we  see  that  the  percent  of  the  TSSQ  explained  by  our  "best"  model,  (TSSQ 
-  RSSQ)/TSSQ,  is  39.3%,  with  a  corresponding  F-statistic,  F^  74  =  3.72  at  a 
significance  level,  p  <  .001.  We  note  that  in  general  the  F-statistic  can 
also  be  used  (albeit  in  an  ad  hoc  manner)  to  measure  the  significance  of 
reductions  in  the  RSSQ  induced  by  adding  parameters  to  a  model,  see  [4J  and 
(5).  Indeed,  we  are  tempted  to  use  this  statistic  to  assess  the  relative 


were  sought  for  the  parameters  0^ 


algorithm  was  then  applied  to  the  "weighting  factors" 
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importance  of  different  mechanisms  of  dispersal  in  our  model  for  the  cabbage 
root  flies.  (After  all,  this  vas  one  of  the  motivations  for  Havkes' 
experiment).  Unfortunately,  the  fit-to-data  of  our  early  models  vas  so  poor 
that  such  a  use  vas  precluded.  In  fact,  ve  observe  that  since  the  calculation 
of  such  statistics  requires  that  TSSQ  >  RSSQ,  F-tests  could  not  be  applied  to 
the  models  of  Tables  4.1,  4.2,  and  4.3.  This  is  because  all  of  the  testing 
through  that  stage  had  resulted  in  values  for  the  RSSQ  that  alvays  exceeded 
198,  vhereas  ve  had  computed  the  TSSQ  to  be  161.  Hovever,  as  ve  have  seen,  ve 
subsequently  identified  a  model,  vhich  provided  not  only  a  statistically 
significant  explanation  of  the  total  sum  of  squares  error,  but  more 
importantly,  a  biologically  meaningful  explanation  of  the  data  in  terms  of 
diffusive,  convective,  and  "grovth/death"  mechanisms  of  dispersal. 

We  close  this  section  vith  a  brief  discussion  of  some  issues  vhich  our 
analysis  raised  concerning  experimental  design  and  the  data  that  results.  A 
characteristic  feature  of  dispersal  experiments  is  the  use  of  a  point  release 
of  the  initial  population.  While  such  a  technique  certainly  facilitates  the 
actual  execution  of  an  experiment,  any  subsequent  mathematical  analysis  must 
deal  vith  the  problem  of  approximating  a  "delta  function"  for  use  as  initial 
data.  But,  more  basic  than  a  purely  mathematical  consideration,  it  has  been 
suggested  that  point  releases  give  rise  to  data  that  over-represent  the  region 
immediately  surrounding  the  release  site  and  under-represent  the  more  distant 
regions  (5],  thus  masking  the  effect  of  any  possible  convective  mechanism. 
Furthermore,  for  the  particular  experiment  ve  considered,  the  use  of  a  point 
release  had  an  additional  impact.  Because  the  release  point  vas  located  at 
the  edge  of  the  array  of  traps,  fully  half  of  the  initial  population  vas  lost 
from  the  model  almost  instantaneously.  As  a  remedy  it  vould  seem  very 
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reasonable  to  try  to  effect  a  distributed  release  of  the  initial  population, 
and  to  do  so  within  the  central  regions  of  the  array. 
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Table  4.1 

Constant  parameter  model 

1)  D  and  a  searched;  =  0.0,  held  fixed 

Estimated:  D  =  349.103 

a  =  51.716 

RSSQ  =  387.98 

2)  Vp  v2,  and  a  searched;  D  =  349.103,  held  fixed 

Estimated:  =  14.915 

v2  =  39.618 

a  =  52.933 

RSSQ  =  308.31 

3)  Vp  v2,  and  a  searched;  D  »  3.49103,  held  fixed* 

Estimated:  =  14.148 

v2  =  31.598 

a  =  59.551 

RSSQ  =  200.63 


[DJ  =  m2/day,  f v1 ]  =  [v2J  =  m/day,  [a]  =  day-1 
N=4 


*Note  the  reduction  of  D  by  two  orders  of  magnitude.  Subsequent  reductions 
did  not  lower  the  RSSQ. 
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Table  4.2 

Time  varying  "grovth/death"  term  added, 

D,  v^,  and  held  fixed  at  final  values  from  Table  4.1 

1)  a^,  k  =  0,  1,  2,  3  searched 

Estimated:  =  59.551 

=*  59.555 
*  21.939 
=  154.200 

RSSQ  =  200.35 

2)  a^,  k  =  0,  1  searched;  =  1.0,  k  =  2,  3  held  fixed 

Estimated:  =  67.773 

«  51.605 

RSSQ  =  197.98 

lajJ  =  day"1 
N=4 


Table  4.3 


Time  varying  convection  terms  added, 

D  and  a(t)  held  fixed  at  final  values  from  Table  4.2 

1)  0^,  Y^,  k  =  2,3  searched;  0^,  y^,  k  =  0,  1  held  fixed  at  final  values 
for  Vj  and  Vj  from  Table  4.1 

Estimated:  =  1007.704 

03  =  282.792 
y2  =  1272.759 
Y3  =  23.053 

RSSQ  =  196.99 

2)  0j*  j  =  0,  1,  2,  3  searched;  Yj  =  2(3^ 

Estimated:  PQ  =  14.148 

01  =  14.148 
02  =  13.501 
03  =  48.060 

RSSQ  =  197.79 

3)  0^.  Yk,  k  =  0,  1  searched;  0^  =  1/20^^  and  Yk  =  l/lvy  k  =  2,3 

Estimated:  0q  =  14.148 

01  =  14.148 
yq  =  31.598 
y:  =  31.598 

RSSQ  =  197.69 


IV  -  IV  =  m/day 
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Table  4.4 


Discontinuous  "growth/death"  term  introduced, 
ot(t)  held  fixed  with  a  discontinuous  profile 


1)  0^  and  k  =  0,  1,  searched;  D 

l/2rlt  k  =  2,3 

Estimated:  0Q  =  14.148 

0J  -  14.148 
yq  =  31.598 
Y1  =  31.598 

RSSQ  =  168.63 


=  3.49  fixed,  ^  =  1/201  and  yk  = 


2)  ^  and  yk,  k  =  2.  3  searched;  D  =  3.49  fixed,  ^  and  Yr,  k  =  0,  1  held 

fixed  at  final  values  for  Vj  and  v2  from  Table  4.1 


Estimated:  02  =  504.519 

03  =  1269.616 


r2  =  14.801 
r3  =  411.220 


RSSQ  =  98.11 


^k  and  rk,  k  =  0,  1  searched;  D  =  3.49  fixed,  ^  and  Yk,  k  =  2,3  held 

_ *1 _  r  * 


fixed  at  values  from  test  2  above 


Estimated: 


0O  =  14.148 
03  =  14.148 
Yq  =  31.598 
YX  =  31.598 


RSSQ  =  98.00 
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4)  D,  0O  and  yQ  searched;  ^  =  eQ,  Yl  =  vQ,  ^  and  yr,  k  =  2,  3  held 

fixed  as  before 

Estimated:  D  =  3.49 

0O  =  14.148 
Y0  =  31.598 

RSSQ  =  98.00 

[D]  =  m2/day,  [0^]  =  [y^]  =  m/dayf  [o^I  =  day-1 
N=4 
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Table  4.5 

"Fine  tuning"  of  convection  terms, 

D  =  3.49  and  a(t)  held  fixed  at  a  discontinuous  profile 

1)  3^  and  yR,  k  =  2,  3  searched;  3^  =  yR  =  31.6,  k  *  0,  1  held  fixed 

Estimated:  32  =  4.07 

33  =  88.5 
Y2  =  75.68 
Y3  =  325.48 

RSSQ  =  105.00 

2)  32  and  33  searched;  all  other  coefficients  held  fixed  at  values  from 
test  1 

Estimated:  32  *  3.424 

33  -  4.143 

RSSQ  =  100.37 

3)  32  and  33  searched;  3^=3^=  14.148  and  remaining  coefficients  held 
fixed  at  values  from  test  1 

Estimated:  32  =  11.01 

33  =  .364 

RSSQ  =  97.73 

[D]  *  m2/day,  [3^]  =  (yr]  =  m/day,  [ojJ  =  day-1 
N=4 
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Table  4.6 

Final  values  of  parameters  for  the  "best"  model 
with  a  continuous  a(t) 

D  =  3.49  m^/day 
3 

v!< t )  ^(t),  where 

k=0 

00  =  14.15  m/day 
0^  =  14.15  m/day 
@2  =  11.01  m/day 
0^  =  *364  m/day 
3 

v2^t)  rklk(t),  where 

k=0  K 

Yq  =  31.6  m/day 
Tj  =  31.6  m/day 
Y2  =  75.68  m/day 
Y^  =  325.48  m/day 
11 

«(t)  =E  where 

k=0 

«q  *  otj  and  0^  a  a  a  . . .  a  ^  and 
0^  =  67.773  day"1 
«2  =  67.773  day'1 
=  61.804  day-1 
0^  =  .003  daj"1 

RSSQ  =  104.71 


Table  4.7 


Spatial  dependence  considered  in  the  convection  terms, 
D  and  a(t)  held  fixed  as  in  Table  4.6 


vx(x) 


E  *vMx> 

k=0  X  * 


v2(y) 


'EVk<y) 

k=0 


Estimated: 


<0  - 

11.851 

Estimated:  Hq 

ss 

10.685 

h  - 

11.871 

* 

19.114 

<2  - 

16.222 

n2 

X 

44.514 

16.152 

n3 

s 

37.451 

II 

.673 

*4 

X 

107.270 

^5 

.885 

n5 

x 

88.090 

^6 

34.978 

*6 

x 

105.171 

<7  " 

35.306 

*7 

x 

65.052 

RSSQ  =  92.42 

iCkJ  =  l\i  =  m/day 
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Table  4.8 

Temporal  dependence  in  convection  terms 
re-estimated  based  on  results  from  Table  4.7 


3 

v1(t,x)  =  (gZ  ^l^Ol-v^x) 

Estimated:  00  =  1.000 

01  =  .998 

02  =  2.425 
03  =  2.519 


RSSQ  =  89.68 


3 

v?(t,y)  =  ( E 

k=0 


Yk1k(t»,v2(y> 


Estimated:  Tq  =  1.001 

YX  =  1.000 

V2  =  1.222 

r3  =  1.455 


13k)  »  l\J  =  m/day 


N=4 
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First  census: 


-.98  -.51  -1.09  1.22 

*  *  *  * 


Figure  4.2 

Residual  error  at  each  grid  point 
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Second  census: 


-.09  .02  .18  .37  -1.46 

*  *  *  *  x  * 


-  y  - 

.7  -2.13 

*  * 


.09 

* 


-.09  .02  .18  -.63  -.46  -.3  -1.13  1.09 

k  k  k  k  k  k 


-.08  .01  .17  .34  .5  -.35  -.2  .01 

★  k  k  k  k  k  k  k 


fallow 


1.22 

-.13 

.51 

.25  .08 

★ 

★ 

★ 

★  ★ 

-.81 

.81 

.49 

.25  .07 

* 

★ 

k 

★  ★ 

2.92 

-.28 

-.56 

.22  -.93 

★ 

★ 

k 

k  k 

.76 

.56 

.32 

.16  -.95 

* 

★ 

k 

★  ★ 

crop 

.38 

.79 

.25 

-1.43  .5 

* 

* 

* 

*  * 

-1.62 

-2.21 

.25 

.57  -1.49 

* 

* 

k 

*  * 

-1.72 

-1.35 

1.08 

.37  -.69 

★ 

* 

★ 

*  * 

.98 

1.26 

1.6 

-1.18  .78 

* 

* 

★ 

★  * 

crop 

The  cross  marks  the  release  point  and  the  arrow  indicates  the  wind  direction. 
RSSQ  =  97.73 


Introduction  of  Random  Error  into  Data 

We  describe  how  we  introduced  random  error  into  our  data  so  that  the 
estimation  scheme  could  be  tested  in  the  presence  of  "noise".  We  first  used 
the  IMSL  routine  GGNML  to  produce  a  set  of  normal  random  numbers  with  mean  0 
and  variance  1  to  correspond  to  our  set  of  analytically  generated  data.  Using 
these  random  numbers  we  perturbed  the  data  poincs  with  the  requirement  that 
the  errors  remain  less  than  ten  percent  with  95%  probability.  That  is,  we 
treated  each  data  point  as  the  mean  of  a  normal  distribution,  then  adjusted 
the  variance  to  insure  that  1.96  standard  deviations  from  the  mean 
corresponded  to  a  ten  percent  deviation  from  the  true  value. 
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